Mesoscale simulations of surfactant dissolution and mesophase formation 
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The evolution of the contact zone between pure surfactant and solvent has been studied by 
mesoscale simulation. It is found that mesophase formation becomes diffusion controlled and fol- 
lows the equilibrium phase diagram adiabatically almost as soon as individual mesophases can be 
identified, corresponding to times in real systems of order 10 /is. 



When pure surfactant comes into contact with water, 
mesophases appear at the interface. This is an impor- 
tant process not only for the practical use of surfactants, 
but also from the point of view of fundamental surfactant 
phase science. Indeed, contact 'flooding' or penetration 
scan experiments can yield quantitative information on 
mesophases as a function of composition 0. The phe- 
nomenon is invariably diffusion controlled in the sense 
that the widths of the mesophases follow t^/'^ growth 
laws, and adiabatic in the sense that the local composi- 
tion determines the mesophase boundaries according to 
the equilibrium phase diagram [|[ |[ ||] . But penetration 
scan experiments are restricted to observation times of 
minutes to hours: what happens on time scales shorter 
than this? How early does diffusion control set in, and 
how soon can one expect the mesophase boundaries to 
track the local composition adiabatically? Such ques- 
tions are not just of scientific interest since time scales 
of seconds or less are important in modern processing 
and the everyday use of surfactants, for instance deter- 
mining how rapidly washing powder dissolves. Experi- 
mentally, this regime is very difficult to access because 
of the short time scales and the relatively small amounts 
of mesophase involved. To probe these questions there- 
fore, we have therefore undertaken novel mesoscale sim- 
ulations of surfactant dissolution. We find adiabatic dif- 
fusion control is established remarkably rapidly in our 
model, on time scales in which only a few repeat units 
of the growing mesophases have appeared, corresponding 
to times in the real systems of order 10 /iS. 

The model we have used is a 'minimalist' particle- 
based model of a binary surfactant / water mix- 
ture, based on the dissipative particle dynamics (DPD) 
method ||. In DPD, the particles are soft spheres, 
interacting with pairwise soft potentials of the form 
U = ^A{1 — r jr^"^ (r < r^) where r is the particle sepa- 
ration, Tc is the range of the interaction, and A the am- 
plitude. In the model we have three species of particles: 
A, B and C. The A and B particles are bound together 
in pairs as dimers with a fixed separation Td, and rep- 
resent the surfactants. The C particles are monomers 
representing the solvent (water). The different species 



are distinguished by their interaction amplitudes, and 
the trick is to find a set of amplitudes which recover 
suitable phase behaviour. Following earlier work we 
use ^AA = ^BB = ^CC = 25, ^AB = 30, Aac = 0, 
^BC — 50 and = 0.5 which gives phase diagram fea- 
tures lying in a convenient temperature range around 
k^T ~ 1 (we fix units by choosing m — Tc — 1 where 
m is the mass of the particles). 

Fig. 1^ shows the phase diagram as a function of 
dimer concentration and fceT, at an overall density 
p — (2Nab + Nc)/V = 6, where dimer concentration 
is defined to be the mole fraction of particles in dimers: 
c = 2 A'^ab/ (2 A^AB + A^c)- Similar to many real sys- 
tems 0, miccllar (Li), hexagonal (Hi) and lamellar (La) 
phases are found in order of increasing dimer concentra- 
tion, and the interactions arc chosen so that on the pure 
dimer side there is an isotropic fluid (L2) phase. The fact 
that a certain amount of solvent is needed to induce the 
La phase reflects the real behaviour of many nonionic 
surfactants, and is crucial to the contact simulations be- 
low. We did not find any cubic phases, although these are 
not always present in real systems (an extensive study of 
lattice models by Larson j|] suggests that cubic phases 
could be engineered to appear particularly if the DPD 
model is elaborated beyond dimers). Also, we caution 
that k^T in the model is not easily mapped onto a real 
temperature. Despite these deficiencies, it is neverthe- 
less remarkable that such a simple model reproduces the 
main features of the phase diagram common to a large 
number of surfactants. 

To simulate surfactant dissolution and mesophase for- 
mation with this model is now very simple. We take two 
simulation boxes, one containing an equilibrated fluid of 
pure dimers and the other containing equilibrated sol- 
vent particles. We place them next to each other and al- 
low the dimers and solvent particles to interdiffuse. Just 
as in the real systems, mesophases start to appear at 
the interface over time. Here we report results for the 
kBT = 1 isotherm, although the kBT = 0.5 isotherm was 
also studied with similar results. We consider simulation 
boxes of size 10^ x 100, with the concentration gradient 
along the long axis. Fig. ^ shows some representative 
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FIG. 1: (a) Phase diagram for dimer / solvent model in (c, fcaT) plane showing points where we find isotropic, hexagonal and 
lamellar phases (crosses, squares and diamonds respectively). The lines are a guide to the eye. (b)-(e) Isosurfaces (pc ~ 1-3) 
for equilibrium phases along ksT = 1 isotherm, at c = 0.3, 0.6, 0.7 and 0.9 respectively, (f)-(h) Isosurfaces showing snapshots 
of surfactant dissolution simulations, at elapsed times t = 12, 120 and 1200 DPD time units from initial contact. 



simulation snapshots. 

Whilst this is straightforward, there are a couple of 
technical points to be considered. Firstly it is impor- 
tant to adjust the densities in the two boxes so that the 
pressures are equal. From separate simulations, we deter- 
mined the equation of state for the pure components, and 
found that densities p = 6.124 and 5.896 for the solvent 
and dimers respectively give a common pressure p « 100. 
These densities are within « 2% of the density p — 6 used 
to construct the phase diagram in Fig. which obviates 
the need to consider constant pressure simulations. The 
pressure matching is important though because it sup- 
presses sound waves which would spoil the subsequent 
analysis. Secondly, if we were to use conventional peri- 
odic boundary conditions, there would be an unwanted 
second interface in the system. We eliminate this by 
bounding the simulation box in the long direction by 
hard reflecting walls, supplemented by short-range soft 
repulsive potentials of the form U = ^Awaii(l — z/rc)"^ 
{z < Tc), where z is the distance from the hard wall. The 
soft repulsive force suppresses density oscillations which 
would otherwise arise due to the abrupt termination of 
the particle density. We find empirically A„aii = 25 gives 
a smoothly vanishing density. Periodic boundary condi- 
tions are retained in the other two directions. 

Now we turn to an objective analysis of the simula- 
tions to determine the growth laws. This is done by in- 
troducing a local order parameter which can distinguish 
between mesophases as a function of distance normal to 
the original contact plane. Although we examined vari- 
ous possibilities such as the use of Minkowski functionals 
, an approach which works well in practice is moti- 



TABLE I: Relationship between the ordered eigenvalues fii of 
the second moment M of the isosurface normal distribution 
p(n), the nature of p(n), the expected isosurface geometry, 
and the expected mesophase; based on Mardia M. 



Eigenvalues of M 


p(n) 


Isosurface 


Mesophase 


m ^ fj,2 ~ /J-s ~ 1/3 


random 


spheres/blobs 


Li or L2 


/il ^ 0, P2 ~ P3 ~ 1/2 


planar 


cylinders 


Hi 


^1 PS /i2 ~ 0, P3 ~ 1 


axial 


planes 


La 



vated by the pc = 1.3 isosurfaces shown in Fig. |l| The 
different mesophases are clearly distinguishable to the eye 
because the isosurface is predominantly spherical, cylin- 
drical, planar, or fragmented, for the Li, Hi, and 
L2 phases respectively. This geometric insight can be 
made quantitative by constructing the second moment 
of the isosurface normal distribution p(n), namely the 
symmetric tensor M = / nnp{n) dn The local geo- 
metric nature of the isosurface and hence the underlying 
mesophase is reflected in the eigenvalues pi of M as laid 
out in Table ||. If these eigenvalues are ranked in order 
of increasing size, we determined after some trials that a 
shce can be classified as Hi if fj,i < 0.05 and p2,3 > 0.3, 
or La if pi < 0.05 and p2 < 0.15 (note J2i=i t^i — !)■ 

For the present study, we divide the simulation box 
into slices of thickness 0.25 Tc parallel to the original con- 
tact plane and determine the eigenvalues of M for each 
slice. Fig. ^(a) gives a representative example, showing 
that the various mesophases can be clearly distinguished. 
The transition between Hi and Lq. is particularly sharp as 



3 




(a) 



-50 50 100 

Distance, z 




1000 
Time, t 



2000 




(c) 



500 1000 1500 2000 

Time, t 



o 6 



<l 



(d) 



' 1 


1 1 1 1 1 












1,1. 



500 1000 1500 
Time, t 



2000 



FIG. 2: (a) Eigenvalues of local order parameter tensor M at t = 1200 DPD time units after initial contact, (b) Phase boundary 
positions as a function of time, as determined from local order parameter tensor eigenvalues, (c) Local concentrations at phase 
boundaries as a function of time, (d) Squared widths of middle mesophases as a function of time. 



can also be seen in Fig. g(h). The above criteria are used 
to classify each slice and the positions of the boundaries 
between mesophases determined as a function of time. 
There is a short incubation period before the local or- 
der parameter can distinguish the different mesophases, 
but beyond this point, the boundaries can be reliably 
tracked as shown in Fig. ||(b). The key results are con- 
tained in Figs. ||(c) and (d), which show respectively the 
local composition at the mesophase boundary and the 
squared width of the middle two mesophases as a func- 
tion of time. 

Fig. ||(c) shows that the local compositions be- 
come established at constant values almost as soon as 
mesophases can be distinguished. Table || shows that 
the Li-Hi and Hi-Lq boundaries in the kinetic simula- 
tion lie almost exactly at the point expected from the 
equilibrium phase diagram. The La-L2 boundary is dis- 
placed from that found in the equilibrium phase diagram, 
but it is clear from Fig. |^(a) that this boundary is not 
completely sharp. Apart from this boundary therefore, it 
appears that the mesophases grow adiabatically almost 
from the earliest moments that the local order parameter 
can distinguish the growing mesophases. 



TABLE IL Composition at mesophase boundaries from equi- 
librium ksT = 1 isotherm in Fig. ^(a), coinpared to those 
determined from kinetic simulation in Fig. Mc). Units are 
dimer concentration c (the figure in brackets is an estimate of 
the error in the final digit). 



Boundary 


Li-Hi 


Hi-Lq 


La-L2 


Equilibrium 


0.435(5) 


0.605(5) 


0.875(5) 


Kinetic 


0.44(1) 


0.61(1) 


0.84(1) 



Fig. ||(d) shows that both mesophases follow a Az^ oc t 
law, where Az is the mesophase width. This is a classic 
verification of diffusion control. To further examine this, 
we fit the slope to what would be expected if the local 
composition followed a simple diffusion law: 

Az^ = 4L>[erf"^(2c2 - 1) - erf"^(2ci - 1)]^^, (1) 

where D is an effective diffusion coefficient, and ci^2 are 
the fixed compositions at the edges of the mesophase of 
interest. 



The results are shown in Table III, where for 
comparison we have included the diffusion coefficients for 
the two pure components, ie dimers in a pure dimer fiuid 
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TABLE III: Slope of A^^ vs t plot from Fig. |(dl, and effective 
diffusion coefficient D determined from Eq. (hi) for the two 
middle mesophases, compared to self diffusion coefficients for 
solvent particles and dimers in the pure components. Units 
are DPD units. 



Mesophase 


Hi 


La 


solvent 


dimers 


Slope 


1.07(1) 


3.23(2) 






Effective D 


3.0(5) 


3.2(4) 


0.219(3) 


0.092(2) 



and solvent particles in pure solvent, determined from 
separate simulations. The effective diffusion coefficients 
are slightly different in the two phases, similar to the 
finding our previous experiments Q , although the differ- 
ence here is within the statistical errors. Thus we con- 
clude that, whilst the diffusion coefficient may be weakly 
composition- and mesophase-dependent, the whole pro- 
cess is diffusion controlled from the earliest times for 
which mesophases can be distinguished. Intriguingly the 
effective diffusion coefficients in the mesophases are con- 
siderably greater than the self diffusion coefficients in the 
pure components. This suggests the existence of a size- 
able free energy driving force enhancing diffusion and is 
reminiscent of a detailed study on the Hi phase of the 
C12E6 / water system ||T^. Note that the orientation 
of the Hi and phases in the simulation is such that 
mutual diffusion occurs in the 'easy' direction: along the 
hexagonal rods, or parallel to the lamellar layers (we have 
also extracted orientation information from the eigenvec- 
tors of M which will be reported elsewhere) . 

We can use the effective diffusion coefficients to map 
elapsed times in our mesoscale simulation onto real times. 
The corresponding diffusion coefficients in nonionic sur- 
factant systems are typically 2 x lO^^^m^s^^ Q. The 
lamellar repeat spacing L can be used as a common mea- 
sure of length: L k, 2.5 DPD units in the simulation, and 
typically L « 5nm in real systems. By using L^/D to 
scale time, we conclude one DPD time step is approx- 
imately equivalent to 50 ns. Thus the duration of the 
whole simulation (2000 DPD time units) represents more 
than O.I ms of real time, equivalent to 10^ molecular dy- 
namics time steps. Fig. ^ shows that adiabatic diffusion 
controlled dissolution sets in after about 200 DPD time 
units, equivalent to about 10 /is, which is the origin of 
the time scale quoted in the introduction. 

Our simulations, whilst addressing very directly in- 
triguing questions concerning surfactant dissolution, are 
also interesting because they are simulations of phase 
formation kinetics starting from a highly inhomogeneous 
state rather than the usual approach which is to quench a 
homogeneous system. To our knowledge, the only com- 



parable study in the past has been the disappearance 
of an interface in a vapour-liquid or binary liquid mix- 
ture after an 'antiquench' (sudden increase in temper- 
ature) It is indeed striking that we do not find 
any evidence of the delayed appearance of mesophases 
due to metastability, unlike that seen experimentally in 
temperature-jump experiments p^ , nor do we find any 
systematic displacement of the mesophase boundaries by 
the strong concentration gradients present at such early 
times. 

Our mesoscale model is also well suited to other 
fundamental (meso)phase studies such as temperature 
quenches 14 , the effect of shear on phase boundaries, 
the exploration of epitaxial relationships such as between 
the Hi and Lq phases in the present simulations, amongst 
many other possibilities. 
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